Method and device for estimating spatial distribution of inertia of power system based on multi-innovation identification

ABSTRACT

A method and device for estimating spatial distribution of power system inertia based on multi-innovation identification. The method includes: (S1) acquiring time-series data including each node frequency of a power system and a transmission power of a transmission line; decomposing a frequency signal in a series of different spaces to identify a disturbance occurring moment; (S2) analyzing a relationship among node inertia, node frequency and node injection power; constructing an OEMA model to describe an active power-frequency dynamic process of the node after the disturbance occurs; analyzing a coupling relationship between an unknown parameter of the model and the node inertia; (S3) solving the OEMA model based on multi-innovation identification to identity the unknown parameter; and (S4) according to the relationship between a vector of the parameter and the node inertia, calculating each node inertia to estimate the spatial distribution of the power system inertia.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims the benefit of priority from Chinese PatentApplication No. 202011525723.5, filed on Dec. 22, 2020. The content ofthe aforementioned application, including any intervening amendmentsthereto, is incorporated herein by reference in its entirety.

TECHNICAL FIELD

The present application relates to dispatching operation of a powersystem, and more particularly to a method and device for estimatingspatial distribution of inertia of a power system based onmulti-innovation identification.

BACKGROUND

With the development of new energy sources and the construction ofhigh-capacity direct current (DC) transmission projects, the synchronousgenerator unit in the power grid is gradually replaced with thenonsynchronous power supply. Since the direct current and new energy areconnected to the grid through a power electronic converter, theirtransmission powers are decoupled from the grid frequency, so that theyare free of inertia response characteristics. With the integration of alarge number of non-synchronous power supplies such as DC and renewableenergies, the overall inertia level of the system will be greatlyattenuated, which will accelerate the frequency change afterencountering a disturbance and deteriorate the frequency stability ofthe system. Additionally, the distribution difference of the new energyresources breaks the original relatively balanced distribution patternof the inertia resources. The centralized grid-connected region of therenewable energy exhibits an obvious low-inertia characteristic, whichis in a sharp contrast with the high-inertia state of the grid-connectedregion of the synchronous generator unit cluster. Additionally, theaccess of the high-capacity DC transmission line blocks the inertiasupport between the grids in different regions to some extent, so thatthe uneven spatial distribution of the system inertia is aggravated. Thereplacement of the local synchronous power supply with the DC feed-inpower further aggravates the imbalance of the inertia distribution ofthe local grid.

With the access of the non-synchronous power supplies, such ashigh-capacity DC and renewable energy to the grid, the power systemexhibits a remarkable inertia spatial distribution characteristic. Underthe active power disturbance, the non-synchronous power supplycentralized region with weak inertia is prone to frequency instabilitydue to the rapider frequency variation with respect to the high-inertiaregions, which will severely limit the overall loading capability of thesystem. Thus, the estimation of the power system inertia should not onlyfocus on the overall level, but also consider the spatial distributionsituation, such that the adaptability of the access of the local newenergy/direct current and the security and stability of the systemduring the operation can be improved.

At present, the estimation of the power system inertia is merely carriedout at the levels of system and region.

1. Estimation of power system inertia based on equivalent inertia

The overall inertia level in the power system is estimated by readingthe on/off status of the units in the supervisory control and dataacquisition (SCADA) system. Specifically, the inertia time constants ofall online units are respectively multiplied by their capacities, andthe products are added to obtain the total rotational kinetic energy ofthe system to determine the overall inertia of the system.

2. Estimation of power system inertia based on polynomial fitting

The regional power grid is simplified as an equivalent synchronous unitto establish an equivalent frequency response model of the regionalpower grid. A frequency simulation curve is obtained by simulating thedisturbance fault, and then the polynomial fitting is employed tofurther extract the initial rate-of-change of frequency (RoCoF) of thesystem after the occurrence of the disturbance. Finally, the overallinertia level of the regional power grid is obtained based on theinitial RoCoF and the capacity of the disturbance fault.

3. Estimation of power system inertia based on system identification

The method commonly performs equivalent simplification on the regionalpower grid, and the frequency deviation, the increment of the loadpower, the power increment of the section tie line, and the increment ofthe mechanical power of the regional power grid are taken as the input.An identification model for describing the active power-frequencyresponse of the regional power grid is constructed to obtain the inertiaof the regional power grid.

From the above analysis, the existing methods for estimating powersystem inertia mainly focus on the inertia level at the system andregion levels, and have limited estimation precision, failing to reflectthe specific spatial distribution of the power system inertia, and alsofailing to provide effective reference information for the operators toidentify the low-inertia region.

SUMMARY

An object of this application is to provide a method and device forestimating spatial distribution of inertia of a power system based onmulti-innovation identification to solve the technical problem that theexisting methods for estimating power system inertia mainly focus on theoverall inertia level at the system and region levels, and have limitedestimation precision, failing to reflect the specific spatialdistribution of the power system inertia and also failing to provideeffective reference information to identify the low-inertia regionduring the operational dispatch of the power system. The device andmethod provided herein can assist the operators to accurately locate theweak-inertia regions of the power system, and arrive at the deepanalysis of the specific distribution of system inertia.

Technical solutions of this application are described as follows.

In a first aspect, this application provides a method for estimatingspatial distribution of inertia of a power system, comprising:

(S1) acquiring a time-series data, wherein the time-series datacomprises a frequency of each node of the power system and atransmission power of a transmission line, and preprocessing thetime-series data; decomposing a frequency signal in a series ofdifferent spaces to accurately identify a moment when a disturbanceoccurs and trace an impact of the disturbance fast;

(S2) analyzing a relationship among a node inertia, a node frequency anda node injection power; constructing an output error moving average(OEMA) model to describe an active power-frequency dynamic process ofeach node after the disturbance; and analyzing a coupling relationshipbetween an unknown parameter in the OEMA model and the node inertia;

(S3) solving the OEMA model by adopting a multi-innovationidentification, and identifying the unknown parameter of the OEMA model;and

(S4) calculating an inertia of each node by using the unknown parameterof the OEMA model identified in step (S3) according to an analyticalrelationship between a vector of the unknown parameter and the nodeinertia to estimate an overall spatial distribution of the inertia ofthe power system.

The working principle of the disclosure is illustrated as follows.

The existing methods for estimating power system inertia mainly focus onthe overall inertia level at the system and region levels, and havelimited estimation precision, failing to reflect the specific spatialdistribution of the power system inertia, and also failing to provideeffective reference information to identify weak-inertia region duringthe dispatching operation process of the power system.

Thus, the disclosure designs a method for estimating spatialdistribution of power system inertia based on multi-innovationidentification, which includes: data preprocessing and accurateidentification of disturbance moments in a research scene; analysis ofthe relationship among the node inertia, node frequency and nodeinjection power and establishment of the OEMA model for estimating thenode inertia; solving of the unknown parameter in the model based on themulti-innovation identification; establishment of a mathematicalrelationship between the unknown parameter and the node inertia; andcalculation of an inertia of each node to further estimate the spatialdistribution of the overall power system inertia.

The data preprocessing and accurate identification of disturbancemoments in the research scene, the rapid tracing of the disturbance, andthe extraction and preprocessing of the data at the inertia responsestage constitute the data foundation for the estimation of the nodeinertia. The analysis of the relationship among the node inertia, nodefrequency, and node injection power and the establishment of the OEMAmodel contribute to converting the estimation of the node inertia intothe identification of the parameters in the model by describing thedynamic change of the node frequency and the node injection power. Bymeans of the multi-innovation identification, the unknown parameter inthe output error moving average (OEMA) model can be quickly solved. Themathematical relationship between the unknown parameter and the nodeinertia is established to calculate the inertia of each node bydeduction, thereby estimating the spatial distribution of the overallpower system inertia. In this application, the estimation of the spatialdistribution of the inertia is refined to the node level, which canprovide the more effective reference information for the dispatchingpersonnel to accurately locate the weak-inertia regions, so as tostrengthen the operation safety and stability of the system and improveits loading capacity.

In an embodiment, the step (S1) comprises:

(S11) acquiring the time-series data for estimating the spatialdistribution of the inertia of the power system in real time using asynchronized phasor measurement unit (PMU), wherein the time-series datacomprises the frequency of each node, the transmission power of thetransmission line, a node load and an electromagnetic power of agenerator;

(S12) preprocessing the time-series data obtained in step (S11), whereinthe preprocessing of the time-series data comprises: interpolation of amissing data and elimination of a data measurement noise adopting afiltering method;

(S13) normalizing the time-series data preprocessed in step (S12) toform a dataset in a standard form for analyzing the node inertia,thereby estimating the spatial distribution of the inertia of the powersystem;

(S14) decomposing the frequency signal in the series of different spacesbased on wavelet multi-resolution analysis using a Daubechies Wavelet;and extracting a wavelet detail coefficient at each level to determine amodulus maxima spot, thereby detecting a singular point of the frequencysignal, wherein the identified singular point of the frequency signal isan initial time when the disturbance occurs.

In an embodiment, the step (S2) comprises:

(S21) in terms of a mechanism of a frequency response of the powersystem, analyzing a coupling relationship among the node inertia, thenode frequency and the node injection power;

(S22) using Laplace transform to convert a time domain relationshipamong the node inertia, the node frequency and the node injection powerinto a continuous frequency domain expression, wherein the continuousfrequency domain expression is a continuous transfer function comprisingan information of the node inertia; and

(S23) constructing the OEMA model to describe the active power-frequencydynamic change process of each node; converting solving of the nodeinertia into identification of the unknown parameter of the OEMA model;and analyzing the coupling relationship between the unknown parameter ofthe OEMA model and the node inertia; wherein an expression of the OEMAmodel is shown as follows:

$\begin{matrix}{{y(t)} = {{\frac{B(z)}{A(z)}{u(t)}} + {{D(z)}{v(t)}}}} \\{{{x(t)} = {\frac{B(z)}{A(z)}{u(t)}}};}\end{matrix}$

A(z)=1+a₁z⁻¹+ . . . +a_(n) _(a) z^(−n) ^(a) ; B(z)=b₁z⁻¹+ . . . +b_(n)_(b) z^(−n) ^(b) ; D(z)=1+d₁z⁻¹+ . . . +d_(n) _(d) z^(−n) ^(d) ; y(t) isan output sequence of the OEMA model to be identified; u(t) is an inputsequence of the OEMA model to be identified; v(t) is a zero-meanGaussian random white noise sequence; x(t) is an internal variable,which represents an actual undetectable variable of the power system,z⁻¹ is a unit backshift operator, expressed as z⁻¹y(t)=y(t−1); and a_(n)_(a) , b_(n) _(b) , d_(n) _(d) represents parameters to be identified ofthe OEMA model;

In an embodiment, in step (S21), the coupling relationship among thenode inertia, the node frequency and the node injection power isexpressed as follows:

${{2H_{i}\frac{d\Delta{f_{i}(t)}}{dt}} = {\Delta{P_{i}(t)}}};$

which is further converted into:

${H_{i} = \frac{\Delta{P_{i}(t)}}{2\frac{d\Delta{f_{i}(t)}}{dt}}};$

wherein H_(i) is an inertia value of node i, and represents a magnitudeof a hindering effect of the inertia of the power system on a frequencychange of the node i during a ΔP_(i)−dΔf_(i)/dt process; Δf_(i) is anode frequency deviation; and ΔP_(i) is an increment of the nodeinjection power;

in step (S22), the continuous transfer function corresponding to therelationship among the node inertia, the node frequency and the nodeinjection power is expressed as follows:

${{G(s)} = {\frac{\Delta{f_{i}(s)}}{\Delta{P_{i}(s)}} = \frac{1}{2H_{i}s}}};$

wherein Δf_(i)(s) is a frequency domain expression of the node frequencydeviation; ΔP_(i)(s) is a frequency domain expression of the incrementof the node injection power; and s is a Laplace operator.

Thus, a pole-zero gain value of a model of the continuous transferfunction is an initial impulse response value, which equals to areciprocal of twice the node inertia.

In an embodiment, the step (S3) comprises:

(S31) expanding a single innovation describing the activepower-frequency dynamic change process of each node into multipleinnovations to establish a multi-innovation matrix to strengthen anidentification precision of the unknown parameter in the OEMA modelestablished in step (S23), wherein the multi-innovation matrix isexpressed as follows:

${\begin{bmatrix}{X\left( {p,t} \right)} \\{Y\left( {p,t} \right)} \\{V\left( {p,t} \right)} \\{\psi\left( {p,t} \right)}\end{bmatrix} = \begin{bmatrix}{{x(t)},{x\left( {t - 1} \right)},\ldots,{x\left( {t - p + 1} \right)}} \\{{y(t)},{y\left( {t - 1} \right)},\ldots,{y\left( {t - p + 1} \right)}} \\{{v(t)},{v\left( {t - 1} \right)},\ldots,{v\left( {t - p + 1} \right)}} \\{{\varphi(t)},{\varphi\left( {t - 1} \right)},\ldots,{\varphi\left( {t - p + 1} \right)}}\end{bmatrix}};$

wherein 1(t) is an information vector comprising the internal variablex(t) and φ(t)=[−x(t−1), . . . , −x(t−n_(a)), u(t−1), . . . , u(t−n_(b)),v(t−1), . . . , v(t−n_(d))]^(T); X(p,t) is an accumulated internalvariable; Y(p,t) is a first accumulated output vector; V(p,t) is anaccumulated noise variable; ψ(p,t) is a first accumulated informationmatrix; and p is the number of an innovation vector;

(S32) replacing an internal undetectable variable of the OEMA model withan output of an auxiliary model based on an auxiliary modelidentification idea using an actual measured information to furthereliminate an impact of noise contained in an input data of the OEMAmodel; wherein the auxiliary model is expressed as follows:{circumflex over (x)}(t)−φ_(S) ^(T)(t)θ_(S);

wherein,

φ_(S)(t)=[−x(t−1), . . . ,−x(t−n_(a)),u(t−1), . . . ,u(t−n_(b))]^(T),θ_(S)=[a₁,a₂, . . . . ,a_(n) _(a) ,b₁,b₂, . . . . ,b_(n) _(b) ] is avector of the parameter to be identified a_(n) _(a) , b_(n) _(b) ; and

(S33) based on step (S32), setting a criterion function; obtaining anoptimal identification result of a parameter of the OEMA model on abasis of a least squares estimation principle on the premise of allowingthe criterion function to be minimum to obtain a value of all unknownparameters of the OEMA model, wherein the criterion function isexpressed as follows:

${{{J(\theta)}:} = {\sum\limits_{i = 1}^{t}{{{Y\left( {p,i} \right)} - {{\psi^{T}\left( {p,i} \right)}\theta}}}^{2}}};$

wherein, Y(p,i) is a second accumulated output vector; ψ(p,i) is asecond accumulated information matrix; φ=[a₁,a₂, . . . . ,a_(n) _(a),b₁,b₂, . . . . ,b_(n) _(b) , d₁, d₂, . . . d_(nd)] represents a vectorof the parameter to be identified a_(n) _(a) , b_(n) _(b) , d_(n) _(d) ;p is the number of the innovation vector.

In an embodiment, the step (S4) comprises:

(S41) considering what the OEMA model describes is a dynamic behavior ofa discrete system, establishing a mapping relationship between a z-planeand an s-plane based on bilinear transformation (Tustin transform); andaccording to the vector of the parameter of the OEMA model identified instep (S3), obtaining an undetermined parameter of a continuous transferfunction corresponding to the OEMA model;

(S42) performing an inverse Laplace transformation on the transferfunction obtained in step (S41) to obtain an impulse response of acorresponding control system; and taking an initial response value toobtain the reciprocal of twice the node inertia, thereby estimating asize of the node inertia; and

(S43) after solving an inertia value of each node of a power grid,estimating the spatial distribution of the inertia of the power system;and visually displaying the spatial distribution of the inertia of thepower system by combining a geographic information of each node.

In an embodiment, in step (S41), the bilinear transformation isconfigured to realize a mapping from the z-plane to the s-plane, and isexpressed as follows:

${z = {\frac{e^{s{T/2}}}{e^{{- s}{T/2}}} \approx \frac{1 + {s{T/2}}}{1 - {s{T/2}}}}};$

wherein, z is a sampling Laplace operator; s is the Laplace operator;and T is a sampling period.

In an embodiment, a discrete transfer function described by the OEMAmodel is expressed as follows:

${{G^{\prime}(z)} = \frac{B(z)}{A(z)}};$

the discrete transfer function is converted into the continuous transferfunction, expressed as follows:

${{G^{\prime}(s)} = \frac{{b_{n - 1}^{\prime}s^{n - 1}} + \ldots + {b_{1}^{\prime}s^{1}} + b_{0}^{\prime}}{{a_{n}^{\prime}s^{n}} + \ldots + {a_{1}^{\prime}s^{1}} + a_{0}^{\prime}}};$

wherein, b_(n-1)′, . . . , b₀′, a_(n)′, . . . , a₀′ are coefficients ofa model of the continuous transfer function; and s is the Laplaceoperator.

In an embodiment, the frequency signal is decomposed and reconstructedunder the Daubechies Wavelet basis to detect the singular point of thefrequency signal to identify an initial moment of a response of theinertia of the power system, wherein an expression of the waveletmulti-resolution analysis is shown as follows:

${{W_{f}\left( {a,b} \right)} = {\left\langle {f,\psi_{a,b}} \right\rangle = {{❘a❘}^{{- 1}/2}{\int_{- \infty}^{+ \infty}{{f(t)}{\psi\left( \frac{t - b}{a} \right)}{dt}}}}}};$

wherein, ψ_(a,b)(t) is a mother wavelet function, which is theDaubechies Wavelet basis; a is a scale coefficient; b is a translatingcoefficient; and f(t) is an acquired frequency signal.

In a second aspect, this application also provides a device forimplementing the above method, comprising:

an acquiring unit;

a first calculating unit;

a second calculating unit; and

an analyzing unit;

wherein the acquiring unit is configured to acquire a time-series dataand preprocess the time-series data, wherein the time-series datacomprises a frequency of each node of a power system and a transmissionpower of a transmission line; the acquiring unit is also configured todecompose a frequency signal in a series of different spaces to identifya moment when a disturbance occurs, and fast trace an impact of thedisturbance;

the first calculating unit is configured to analyze a relationship amonga node inertia, a node frequency and a node injection power, constructan output error moving average (OEMA) model to describe anactive-power-frequency dynamic process of each node after thedisturbance, and analyze a coupling relationship between an unknownparameter of the OEMA model and the node inertia;

the second calculating unit is configured to solve the OEMA model in thefirst calculating unit by adopting a multi-innovation identification toidentify the unknown parameter of the OEMA model; and

the analyzing unit is configured to calculate an inertia of each nodeusing the unknown parameter of the OEMA model identified by the secondcalculating unit according to an analytical relationship between avector of the unknown parameter and the node inertia to estimate anoverall spatial distribution of the inertia of the power system.

In an embodiment, the acquiring unit is operated through the followingsteps:

acquiring the time-series data for estimating the spatial distributionof the inertia of the power system in real time using a synchronizedphasor measurement unit (PMU), wherein the data comprises the frequencyof each node, the transmission power of the transmission line, a nodeload and an electromagnetic power of a generator;

preprocessing the acquired time-series data, wherein the preprocessingof the acquired time-series data comprises: interpolation of a missingdata and elimination of a data measurement noise by using a filteringmethod;

normalizing the preprocessed time-series data to generate a dataset in astandard form for analyzing the node inertia, thereby estimating thespatial distribution of the inertia of the power system; and

decomposing the frequency signal in the series of different spaces basedon wavelet multi-resolution analysis using a Daubechies Wavelet; andextracting a wavelet detail coefficient at each level to determine amodulus maxima spot, thereby detecting a singular point of the frequencysignal, wherein the identified singular point of the frequency signal isan initial time when the disturbance occurs.

In an embodiment, the first calculating unit is operated through stepsof:

in terms of a mechanism of a frequency response of the power system,analyzing a coupling relationship among the node inertia, the nodefrequency and the node injection power; wherein the couplingrelationship among the node inertia, the node frequency and the nodeinjection power is expressed as follows:

${{2H_{i}\frac{d\Delta{f_{i}(t)}}{dt}} = {\Delta{P_{i}(t)}}};$

which is further converted into:

${H_{i} = \frac{\Delta{P_{i}(t)}}{2\frac{d\Delta{f_{i}(t)}}{dt}}};$

wherein H_(i) is an inertia value of node i, and represents a magnitudeof a hindering effect of the inertia of the power system on a frequencychange of the node i during a ΔP_(i)−dΔf_(i)/dt process; Δf_(i) is anode frequency deviation; and ΔP_(i) is an increment of the nodeinjection power;

using Laplace transform to convert a time domain relationship among thenode inertia, the node frequency and the node injection power into acontinuous frequency domain expression, wherein the continuous frequencydomain expression is a continuous transfer function comprising aninformation of the node inertia, and the continuous transfer function isexpressed as follows:

${{G(s)} = {\frac{\Delta{f_{i}(s)}}{\Delta{P_{i}(s)}} = \frac{1}{2H_{i}s}}};$

wherein Δf_(i)(s) is a frequency domain expression of the node frequencydeviation; ΔP_(i)(s) is a frequency domain expression of the incrementof the node injection power; and s is a Laplace operator;

thus, a pole-zero gain value of a model of the continuous transferfunction is an initial impulse response value, which equals to areciprocal of twice the node inertia;

constructing the OEMA model to describe the active power-frequencydynamic change process of each node; converting solving of the nodeinertia into identification of the unknown parameter of the OEMA model;and analyzing the coupling relationship between the unknown parameter ofthe OEMA model and the node inertia; wherein an expression of the OEMAmodel is shown as follows:

$\begin{matrix}{{y(t)} = {{\frac{B(z)}{A(z)}{u(t)}} + {{D(z)}{v(t)}}}} \\{{{x(t)} = {\frac{B(z)}{A(z)}{u(t)}}};}\end{matrix}$

wherein A(z)=1+a₁z⁻¹+ . . . +a_(n) _(a) z^(−n) ^(a) ; B(z)=b₁z⁻¹+ . . .+b_(n) _(b) z^(−n) ^(b) ; D(z)=1+d₁z⁻¹+ . . . +d_(n) _(d) z^(−n) ^(d) ;y(t) is an output sequence of the OEMA model to be identified; u(t) isan input sequence of the OEMA model to be identified; v(t) is azero-mean Gaussian random white noise sequence; x(t) is an internalvariable, which represents an actual undetectable variable of the powersystem; z⁻¹ is a unit backshift operator, expressed as z⁻y(t)=y(t−1);and a_(n) _(a) , b_(n) _(b) , d_(n) _(d) represent parameters to beidentified of the OEMA model.

In an embodiment, the second calculating unit is operated through stepsof:

expanding a single innovation describing the active power-frequencydynamic change process of each node into multiple innovations togenerate a multi-innovation matrix to strengthen an identificationprecision of the unknown parameter in the OEMA model established by thefirst calculating unit, wherein the multi-innovation matrix is expressedas follows:

${\begin{bmatrix}{X\left( {p,t} \right)} \\{Y\left( {p,t} \right)} \\{V\left( {p,t} \right)} \\{\psi\left( {p,t} \right)}\end{bmatrix} = \begin{bmatrix}{{x(t)},{x\left( {t - 1} \right)},\ldots,{x\left( {t - p + 1} \right)}} \\{{y(t)},{y\left( {t - 1} \right)},\ldots,{y\left( {t - p + 1} \right)}} \\{{v(t)},{v\left( {t - 1} \right)},\ldots,{v\left( {t - p + 1} \right)}} \\{{\varphi(t)},{\varphi\left( {t - 1} \right)},\ldots,{\varphi\left( {t - p + 1} \right)}}\end{bmatrix}};$

wherein φ(t) is an information vector comprising the internal variablex(t) and φ(t)=[−x(t−1), . . . , −x(t−n_(a)), u(t−1), . . . , u(t−n_(b)),v(t−1), . . . , v(t−n_(d))]^(T); X(p,t) is an accumulated internalvariable; Y(p,t) is a first accumulated output variable; V(p,t) is anaccumulated noise variable; ψ(p,t) is a first accumulated informationmatrix; and p is the number of an innovation vector;

replacing an internal undetectable variable of the OEMA model with anoutput of an auxiliary model based on an auxiliary model identificationidea using an actual measured information to further eliminate an impactof noise contained in an input data of the OEMA model; wherein theauxiliary model is expressed as follows:{circumflex over (x)}(t)−φ_(S) ^(T)(t)θ_(S);

wherein, φ_(S)(t)=[−x(t−1), . . . ,−x(t−n _(a)),u(t−1), . . . ,u(t−n_(b))]^(T), θ_(S)=[a₁,a₂, . . . . ,a_(n) _(a) ,b₁,b₂, . . . . ,b_(n)_(b) ] is a vector of the parameter to be identified a_(n) _(a) , b_(n)_(b) ; and

setting a criterion function; obtaining an optimal identification resultof the parameter of the OEMA model on a basis of a least squaresestimation principle on the premise of allowing the criterion functionto be minimum, to obtain a value of all unknown parameters of the OEMAmodel, wherein the criterion function is expressed as follows:

${{{J(\theta)}:} = {\sum\limits_{i = 1}^{t}{{{Y\left( {p,i} \right)} - {{\psi^{T}\left( {p,i} \right)}\theta}}}^{2}}};$

wherein, Y(p,i) is a second accumulated output vector; ψ(p,i) is asecond accumulated information matrix; φ=[a₁,a₂, . . . . ,a_(n) _(a),b₁,b₂, . . . . ,b_(n) _(b) , d₁, d₂, . . . , d_(nd)] represents avector of the parameter to be identified a_(n) _(a) , b_(n) _(b) , d_(n)_(d) ; p is the number of the innovation vector.

In an embodiment, the analyzing unit is operated through steps of:

considering what the OEMA model built by the first calculating unitdescribes is a dynamic behavior of a discrete system, establishing amapping relationship between a z-plane and an s-plane based on bilineartransformation (Tustin method); and according to the vector of theparameter of the OEMA model identified by the second calculating unit,obtaining an undetermined parameter of a continuous transfer functioncorresponding to the OEMA model;

performing an inverse Laplace transformation on the obtained transferfunction to obtain an impulse response of a corresponding controlsystem; and taking an initial response value to obtain the reciprocal oftwice the nodal inertia, thereby estimating a size of the node inertia;and

after solving an inertia value of each node of a power grid, estimatingthe spatial distribution of the inertia of the power system; andvisually displaying the spatial distribution of the inertia of the powersystem by combining a geographic information of each node.

In an embodiment, the bilinear transformation is configured to realize amapping from the z-plane to the s-plane, and is expressed as follows:

${z = {\frac{e^{s{T/2}}}{e^{{- s}{T/2}}} \approx \frac{1 + {s{T/2}}}{1 - {s{T/2}}}}};$

wherein, z is a sampling Laplace operator; s is the Laplace operator;and T is a sampling period.

In an embodiment, a discrete transfer function described by the OEMAmodel is expressed as follows:

${{G^{\prime}(z)} = \frac{B(z)}{A(z)}};$

the discrete transfer function is converted into the continuous transferfunction, expressed as follows:

${{G^{\prime}(s)} = \frac{{b_{n - 1}^{\prime}s^{n - 1}} + \ldots + {b_{1}^{\prime}s^{1}} + b_{0}^{\prime}}{{a_{n}^{\prime}s^{n}} + \ldots + {a_{1}^{\prime}s^{1}} + a_{0}^{\prime}}};$

wherein, b_(n-1)′, . . . , b₀′, a_(n)′, . . . , a₀′ are coefficients ofa model of the continuous transfer function; and s is the Laplaceoperator.

In an embodiment, the frequency signal is decomposed and reconstructedunder the Daubechies Wavelet to detect the singular point of thefrequency signal to identify an initial moment of a response of theinertia of the power system, wherein an expression of the waveletmulti-resolution analysis is shown as follows:

${{W_{f}\left( {a,b} \right)} = {\left\langle {f,\psi_{a,b}} \right\rangle = {{❘a❘}^{{- 1}/2}{\int_{- \infty}^{+ \infty}{{f(t)}{\psi\left( \frac{t - b}{a} \right)}{dt}}}}}};$

wherein, ψ_(a,b)(t) is a mother wavelet function, and herein is theDaubechies Wavelet; a is a scale coefficient; b is a translatingcoefficient; and f(t) is an acquired frequency signal.

Compared with the prior art, this disclosure has the followingbeneficial effects.

1. The method provided herein can offer a precise estimation for thespatial distribution of the power system inertia. Compared with theconventional inertia estimation methods, the method provided hereinestablishes an analytical relationship among the node inertia, the nodefrequency and the node injection power, and estimates the node inertiabased on the output error moving average (OEMA) model, thereby refiningthe estimation of the inertia spatial distribution to a node level. Byestimating the inertia of each node in the system, the spatialdistribution of the power system inertia can be obtained, facilitatinghelping the dispatching personnel accurately locate the weak-inertiaregions, and arriving at the deep analysis of the distribution of systeminertia.

2. The model constructed herein for estimating the spatial distributionof the inertia is of strong robustness. The OEMA model constructed fordescribing the active power-frequency dynamic change process of the nodeis a high-order model, in which the time-series data information, suchas the node frequency, is fully considered. Compared with a low-ordermodel merely using the node data at a single moment for estimating thenode inertia, the high-order model can effectively overcome theinterference of a large number of random factors such as datameasurement noise, and avoid being excessively affected by the accuracyof data measurement at a single moment, allowing for a strongrobustness.

3. The model constructed herein for estimating the spatial distributionof the inertia can be fast and effectively solved. Considering that theconstructed model is a high-order model, the single innovation isexpanded to the multiple innovations to match the model, and by means ofthe multi-innovation identification, the identification accuracy isenhanced, so that the node inertia can be fast and effectively estimatedin a short time, analyzing the spatial distribution of the inertia ofthe system.

4. The model for estimating the inertia spatial distribution hasdesirable adaptability for the disturbance fault. Since the constructedOEMA model comprehensively considers the coupling relationship betweenthe node inertia, the node frequency and the node injection power, anduses the time-series information of the node data, it can effectivelyovercome the impact of the disturbance fault on the estimation result ofthe inertia spatial distribution, achieving the accurate estimation ofthe spatial distribution of the system inertia under variousdisturbances (varying in location and intensity).

BRIEF DESCRIPTION OF THE DRAWINGS

The patent or application file contains at least one drawing executed incolor. Copies of this patent or patent application publication withcolor drawing(s) will be provided by the Office upon request and paymentof the necessary fee.

The drawings described herein are intended to facilitate theunderstanding of the application, which are merely illustrative and arenot intended to limit the disclosure.

FIG. 1 is a flowchart of a method for estimating spatial distribution ofan inertia of a power system based on multi-innovation identificationaccording to an embodiment of this disclosure.

FIG. 2 is a topological structure diagram of an IEEE-39 bus systemaccording to an embodiment of this disclosure.

FIG. 3 schematically illustrates estimation results of a node inertiaunder disturbances varying in location and intensity according to anembodiment of this disclosure.

FIG. 4 is a visual display diagram of the inertia spatial distributionobtained by the method of this disclosure.

FIG. 5 is a visual display diagram of the inertia level of the regionalpower grid obtained based on a conventional method.

DETAILED DESCRIPTION OF EMBODIMENTS

The disclosure will be described in detail below with reference to theembodiments and accompanying drawings to make the objectives, technicalsolutions, and advantages of the disclosure clearer. It should beunderstood that these embodiments are only illustrative of thedisclosure, and are not intended to limit the disclosure.

A large number of specific details are described below to provide athorough understanding of this disclosure. However, it is obvious forthose skilled in the art that the present application can be implementedwithout some of these specific details. In other embodiments, thestructures, electric circuits, materials and methods commonly known inthe art are not specifically described to avoid obscuring thisdisclosure.

Throughout the description, references to “one embodiment”, “anembodiment”, “one example”, or “an example” mean that a specificfeature, structure, or characteristic described in conjunction with theembodiment or example is included in at least one embodiment of thisdisclosure. Thus, as used herein, the phrase “one embodiment”, “anembodiment”, “one example” or “an example” does not necessarily refer tothe same embodiment or example. Additionally, specific features,structures, or characteristics can be combined in any appropriatecombination and/or sub-combination in one or more embodiments orexamples. Apart from that, it should be understood by those of ordinaryskill in the art that the drawings provided herein are onlyillustrative, and the drawings are not necessarily drawn to scale. Theterm “and/or” used herein is intended to include any and allcombinations of one or more of the listed items.

As used herein, it should be understood that the orientation or positionrelationship indicated by terms “front”, “rear”, “left”, “right”,“upper”, “lower”, “vertical”, “horizontal”, “high”, “low”, “inner” and“outer” is based on the orientation or position relationship shown inthe accompanying drawings, and is only for the convenience of describingthis disclosure and simplifying the description, rather than indicatingor implying that the referred device or element must have a specificorientation or be constructed or operated in a specific orientation.Therefore, these terms cannot be understood as a limitation for theprotection scope of this disclosure.

Embodiment 1

Referring to an embodiment shown in FIGS. 1-5, a method for estimatingspatial distribution of an inertia of a power system is illustrated,which includes the following steps.

(S1) A time-series data is acquired and preprocessed, where thetime-series data includes a frequency of each node of the power systemand a transmission power of a transmission line. A frequency signal isdecomposed in a series of different spaces to accurately identify amoment when a disturbance occurs and rapidly trace an impact of thedisturbance.

(S2) A relationship among a node inertia, a node frequency, a nodeinjection power is analyzed. An output error moving average (OEMA) modelis constructed to describe an active power-frequency dynamic process ofeach node after the disturbance, and a coupling relationship between anunknown parameter in the OEMA model and the node inertia is analyzed.

(S3) The OEMA model is solved by adopting a multi-innovationidentification to identify the unknown parameter of the OEMA model.

(S4) An inertia of each node is calculated by using the unknownparameter of the OEMA model identified in step (S3) according to ananalytical relationship between a vector of the unknown parameter andthe node inertia to estimate a spatial distribution of the inertia ofthe power system.

In this embodiment, the acquisition and preprocess of the data are asfollows.

A data for estimating the spatial distribution of the inertia of thepower system in real time is acquired by using a synchronized phasormeasurement unit (PMU), in which the data includes the frequency of eachnode, the transmission power of the transmission line, a node load andan electromagnetic power of a generator. The obtained data ispreprocessed, in which the preprocessing of the data includes:interpolation of a missing data and elimination of a data measurementnoise adopting a filtering method. The preprocessed time-series data isnormalized to form a dataset in a standard form for analyzing the nodeinertia, thereby estimating the spatial distribution of the inertia ofthe power system.

The identification of the moment when the disturbance occurs isdescribed as follows.

The frequency signal is decomposed in the series of different spacesbased on wavelet multi-resolution analysis using a Daubechies Wavelet.And a wavelet detail coefficient at each level is extracted to determinea modulus maxima spot, thereby detecting a singular point of thefrequency signal, in which the identified singular point of thefrequency signal is an initial time when the disturbance occurs. Anexpression of a decomposition of the frequency signal based on theDaubechies Wavelet basis is shown as follows:

${{W_{f}\left( {a,b} \right)} = {\left\langle {f,\psi_{a,b}} \right\rangle = {{❘a❘}^{{- 1}/2}{\int_{- \infty}^{+ \infty}{{f(t)}{\psi\left( \frac{t - b}{a} \right)}{dt}}}}}};$

where ψ_(a,b)(t) is a mother wavelet function, which is the DaubechiesWavelet basis; a is a scale coefficient; b is a translating coefficient;and f(t) is an acquired frequency signal.

The construction of the model of estimation of the node inertia is asfollows.

The node inertia is used to characterize the intensity of the inertia atdifferent electrical positions of the power system by estimating eachnode inertia value of the power system, so as to assess the specificsituation of the spatial distribution of the power system. The nodeinertia is an inherent property of system inertia that hinders the nodefrequency change during the process of the energy fluctuation in thepower system, which measures the ability of the node to resist activepower disturbance under the combined action of all inertia resources. Interms of a mechanism of a frequency response of the power system, thecoupling relationship among the node inertia, the node frequency and thenode injection power at the stage of the inertia response after thedisturbance fault is expressed as follows:

$\begin{matrix}{{{2H_{i}\frac{d\Delta{f_{i}(t)}}{dt}} = {\Delta{P_{i}(t)}}};} & (2.1)\end{matrix}$

which is further converted into:

$\begin{matrix}{{H_{i} = \frac{\Delta{P_{i}(t)}}{2\frac{d\Delta{f_{i}(t)}}{dt}}};} & (2.2)\end{matrix}$

where H_(i) is an inertia value of node i, and represents a magnitude ofa hindering impact of the inertia of the power system on a frequencychange of the node i during a transient process; Δf_(i) is a nodefrequency deviation; ΔP_(i) is an increment of the node injection power.

The continuous transfer function corresponding to the relationship amongthe node inertia, the node frequency and the node injection power inequation (2.1) is expressed as follows:

$\begin{matrix}{{G(s)} = {\frac{\Delta{f_{i}(s)}}{\Delta{P_{i}(s)}} = {\frac{1}{2H_{i}s}.}}} & (3)\end{matrix}$

It can be seen from equation (3) that a pole-zero gain value of a modelof the continuous transfer function is an initial impulse responsevalue, which equals to a reciprocal of twice the node inertia H_(i)mathematically.

As equation (2.1) shows, the node inertia information reflected duringthe process of the node active power-frequency dynamic change, andequation (3) shows that if the pole-zero gain value of G(s) is obtained,the node inertia is estimated. In order to avoid being excessivelyaffected by the accuracy of data measurement at a single moment, inwhich an expression of the OEMA model describing the process is shown asfollows:

$\begin{matrix}{{y(t)} = {{\frac{B(z)}{A(z)}{u(t)}} + {{D(z)}{v(t)}}}} \\{{{x(t)} = {\frac{B(z)}{A(z)}{u(t)}}};}\end{matrix}$

in which, A(z)=1+a₁z⁻¹+ . . . +a_(n) _(a) z^(−n) ^(a) , B(z)=b₁z⁻¹+ . .. +b_(n) _(b) z^(−n) ^(b) , D(z)=1+d₁z⁻¹+ . . . +d_(n) _(d) z^(−n) ^(d), y(t) is an output sequence of a model to be identified (OEMA model);u(t) is an input sequence of the model to be identified (OEMA model);v(t) is a zero-mean Gaussian random white noise sequence; x(t) is aninternal variable, which represents an actual undetectable variable ofthe power system, z⁻¹ is a unit backshift operator expressed asz⁻¹y(t)=y(t−1); and a_(n) _(a) , b_(n) _(b) , d_(n) _(d) representsparameters to be identified of the OEMA model.

The identification based on the multi-innovation identification of theunknown parameter in the OEMA model is as follows.

From the above, the unknown parameter in the OEMA model is firstlyidentified to assess the node inertia. The innovation is a usefulinformation for improving the estimation accuracy of the parameters orthe state estimation accuracy of the state. The multi-innovationidentification expands the original single innovation to multipleinnovations. With the use of the variable length of the innovation, thenode data measured by the phasor measurement unit (PMU) is fullyutilized to improve the utilization efficiency of the data, andstrengthen the identification precision of the system. In themulti-innovation identification, p (p is the number of innovations)groups of the data in the data windows from (t−p+1) to t, an expandedinternal variable x(t), an output vector y(t), an noise vector v(t), aninformation vector φ(t) are considered to establish a multi-innovationmatrix, which is shown as follows:

$\begin{matrix}{\begin{bmatrix}{X\left( {p,t} \right)} \\{Y\left( {p,t} \right)} \\{V\left( {p,t} \right)} \\{\psi\left( {p,t} \right)}\end{bmatrix} = \begin{bmatrix}{{x(t)},{x\left( {t - 1} \right)},\ldots,{x\left( {t - p + 1} \right)}} \\{{y(t)},{y\left( {t - 1} \right)},\ldots,{y\left( {t - p + 1} \right)}} \\{{v(t)},{v\left( {t - 1} \right)},\ldots,{v\left( {t - p + 1} \right)}} \\{{\varphi(t)},{\varphi\left( {t - 1} \right)},\ldots,{\varphi\left( {t - p + 1} \right)}}\end{bmatrix}} & (5)\end{matrix}$

in which, φ(t) is an information vector including the internal variablex(t), and φ(t)=[−x(t−1), . . . , −x(t−n_(a)), u(t−1), . . . ,u(t−n_(b)), v(t−1), . . . , v(t−n_(d))]^(T); X(p,t) is an accumulatedinternal variable, Y(p,t) is a first accumulated output vector, v(p,t)is an accumulated noise variable, ψ(p,t) is a first accumulatedinformation matrix.

The frequency data actually acquired by the phasor measurement unit(PMU) usually includes the random noise, which is difficult to fullyreflect the actual frequency dynamic process of the node. Therefore, theauxiliary model identification idea is adopted. With the assist of themeasurable information of the system, the output of the auxiliary model{circumflex over (x)}(t) replaces the undetectable variable x(t) in thesystem to further strengthen the identification precision of the unknownparameters in the OEMA model. The auxiliary model is expressed asfollows:{circumflex over (x)}(t)−φ_(S) ^(T)(t)θ_(S)  (6);

wherein, φ_(S)(t)=[−x(t−1), . . . ,−x(t−n_(a)),u(t−1), . . .,u(t−n_(b))]^(T), and θ_(S)=[a₁,a₂, . . . . ,a_(n) _(a) ,b₁,b₂, . . . .,b_(n) _(b) ] is a vector of the parameter to be identified.

Auxiliary model-based multi-innovation extended least-squares algorithm(AM-MI-ELS algorithm, hereinafter referred to as “multi-innovationidentification”) is adopted to take the node injection power as aninput, and the node frequency as an output, so as to update theestimation results of the unknown parameters constantly and fast, andeach fast updated result is expressed as follows:

$\begin{matrix}{{\overset{\hat{}}{\theta} = {\left\lbrack {\sum\limits_{i = 1}^{t}{{\psi\left( {p,i} \right)}{\psi^{T}\left( {p,i} \right)}}} \right\rbrack^{- 1}\left\lbrack {\sum\limits_{i = 1}^{t}{{\psi\left( {p,i} \right)}{Y^{T}\left( {p,i} \right)}}} \right\rbrack}};} & (7)\end{matrix}$

the prediction error is made as less as possible, and the settingcriterion function is as follows:

$\begin{matrix}{{{{J(\theta)}:} = {\sum\limits_{i = 1}^{t}{{{Y\left( {p,i} \right)} - {{\psi^{T}\left( {p,i} \right)}\theta}}}^{2}}};} & (8)\end{matrix}$

based on the optimal estimation of the least squares principle, when theidentification result is updated till the criterion function is minimum,the unknown parameters in the OEMA model is identified.

The estimation of the node inertia is as follows:

The OEMA model is a discrete model, and the discrete transfer functioncorresponding to it is as follows:

$\begin{matrix}{{{G^{\prime}(z)} = \frac{B(z)}{A(z)}};} & (9)\end{matrix}$

a mapping relationship between a z-plane and an s-plane is constructedto convert the discrete transfer function into the continuous transferfunction based on bilinear transform (Tustin's transform), and isexpressed as follows:

$\begin{matrix}{{z = {\frac{e^{s{T/2}}}{e^{{- s}{T/2}}} \approx \frac{1 + {s{T/2}}}{1 - {s{T/2}}}}};} & (10)\end{matrix}$

the continuous transfer function after the transform is expressed asfollows:

$\begin{matrix}{{G^{\prime}(s)} = \frac{{b_{n - 1}^{\prime}s^{n - 1}} + \ldots + {b_{1}^{\prime}s^{1}} + b_{0}^{\prime}}{{a_{n}^{\prime}s^{n}} + \ldots + {a_{1}^{\prime}s^{1}} + a_{0}^{\prime}}} & (11)\end{matrix}$

-   where b_(n-1)′, . . . , b₀′, a_(n)′, . . . , a₀′ are undetermined    coefficients of a model of the continuous transfer function; and s    is the Laplace operator.

The discrete transfer function G′(z) corresponding to the OEMA model isperformed with the Laplace transform to obtain a continuous transferfunction G′(s) and G′(s) is performed with the inverse Laplace transformto obtain an impulse response of a corresponding control system. Theinitial response value is taken to obtain a reciprocal of twice the nodeinertia, thereby estimating the magnitude of the node inertia.

After an inertia value of each node of the power grid is solved, thespatial distribution of the inertia of the power system is estimated andexpressed as follows:H={H _(i) |i=1,2, . . . ,n}   (12);

in which, H is a collection of all node inertia of the power grid,reflecting the specific situation of the spatial distribution of theinertia of the power system, n is the number of the nodes of the powersystem.

On the basis of equation (12), the information of the node inertia sizeand the node geographic location obtained by fusion calculation arevisualized to realize the visual display of the spatial distribution ofthe inertia of the power system.

Based on the above, the method for estimation of spatial distribution ofthe inertia of the power system based on multi-innovation identificationis applied to the New England IEEE 39-Bus System to assess the spatialdistribution of the inertia of the power system. The calculation resultof the frequency stability simulates the synchronization measurementresult by the phasor measurement unit (PMU) in the actual system. InFIG. 2, the topological structure diagram of the IEEE 39 bus system isshown, and 10 generators are provided herein, with the total load of thesystem of 6150 MW and the rated frequency of 50 Hz. The model parametersused by the generator are shown in Table 1.

TABLE 1 Parameters of the generator models Generator H x′_(d) x′_(q)x_(d) x_(q) T′_(d0) T′_(q0) x_(l) 1 5.95 0.060 0.080 0.200 0.19 7.000.70 0.030 2 2.58 0.697 1.700 2.950 2.82 6.56 1.50 0.350 3 3.90 0.5310.876 2.495 2.37 5.70 1.50 0.304 4 3.20 0.436 1.660 2.620 2.58 5.69 1.500.295 5 3.62 1.320 1.660 6.700 6.20 5.40 0.44 0.540 6 3.79 0.500 0.8142.540 2.41 7.30 0.40 0.224 7 3.34 0.490 1.860 2.950 2.92 5.66 1.50 0.3228 3.18 0.570 0.911 2.900 2.80 6.70 0.41 0.280 9 2.94 0.570 0.587 2.1062.05 4.79 1.96 0.298 10 3.54 0.310 0.080 1.000 0.69 10.20 0 0.125

During the simulation analysis, to analyze the adaptability of theproposed methods for estimation of the spatial distribution of powersystem inertia to different disturbance faults varying in intensity, thefault scenarios varying in intensity are set as follows.

In a large disturbance scenario, the specific setting of thelarge-disturbance fault is as follows. At 10 s, a total load of 615 MW(approximately 10% of the total load) is removed at bus #20, to verifythe effectiveness of this method in estimating the spatial distributionof the inertia under the condition of large disturbance.

In a small disturbance scenario, the specific setting of thesmall-disturbance fault is as follows. At nodes of #4, #7, #15, #20,#23, #27, the small load disturbance of 100 MW (approximately 1.63% ofthe total load) is added to simulate the small disturbance fault such asthe load switching during the operation of the stable state of thesystem, so as to verify the effectiveness of this method in estimatingthe spatial distribution of the inertia under the conditions of smalldisturbance or quasi-stable-state.

The general nodes are considered to have no accurate inertia value forcomparison and verification, so the active power-frequency data of thegenerator node after the fault is selected to estimate the inertia ofeach synchronous unit. The results are shown in Table 2. The inertiatime constant of the generator is based on its rated capacity.

TABLE 2 Estimation results of inertia time constant of generators Gen-Actual Large disturbance fault Small disturbance fault erator inertia ofEstimated Erroll Estimated Erroll number generator/s value/s % value/s %G1 5.95 6.10 2.61 5.82 2.24 G2 2.58 2.68 4.10 2.47 3.97 G3 3.90 3.971.90 3.75 3.76 G4 3.20 3.31 3.23 3.06 4.40 G5 3.62 3.70 2.08 3.55 1.98G6 3.79 3.72 1.90 3.62 4.49 G7 3.34 3.28 1.84 3.40 1.74 G8 3.18 3.034.97 3.14 1.30 G9 2.94 2.98 1.08 3.11 5.62 G10 3.54 3.60 1.74 3.72 5.02

As shown in Table 2, the multi-innovation identification is adopted tocompare the estimated inertia value of each unit with the actual inertiavalue. Under the condition of large disturbance fault, the estimatederror of the inertia of each unit is within 5%, the minimum error is1.08%, and the maximum error is 4.97%. Under condition of smalldisturbance fault, the estimated error of the inertia of each unit iswithin 6%, the minimum error is 1.3%, and the maximum error is 5.62%. Itillustrates that the selection of the data obtained after thedisturbances varying in intensity for identifying the inertia of thegenerator has almost no impact on the estimation result, which verifiesthat the method for estimation of the spatial distribution of inertiabased on the multi-innovation identification herein accurately analyzesthe spatial distribution of the power system inertia under disturbancesvarying in intensity.

In this embodiment, the disturbance faults varying in intensity andlocation are set to verify the adaptability of the method for estimationof the spatial distribution of the inertia to disturbances varying inlocation. During the simulation analysis, the settings of the fault areas follows.

Scenario 1: at 10 s, the load at node #3 is cut off to 200 MW(approximately 3.25% of the total load).

Scenario 2: at 10 s, the load at node #8 is cut off to 200 MW(approximately 3.25% of the total load).

Scenario 3: at 10 s, the load at node #15 is cut off to 200 MW(approximately 3.25% of the total load).

Scenario 4: at nodes of #4, #7, #15, #20, #23, #27, a small loaddisturbance of 100 MW (about 1.63% of the total load) is added tosimulate the fluctuation of the overall active power in the systemduring the normal operation.

In FIG. 3, the estimation results of a node inertia under disturbancesvarying in location are given. It can be seen from FIG. 3 that theestimated values of the node inertia under disturbances varying inlocation and intensity based on the method provided herein are veryclose. Among the 39 nodes of the overall system, the deviation betweenthe maximum and minimum values of the estimation results of the inertiais within 0.4 s, and the estimated result of the inertia is of nooutlier. It illustrates that the model of the estimation of the spatialdistribution of the inertia herein is of good adaptability to differentdisturbances and strong robustness.

The IEEE 39 bus system topology is mapped to the two-dimensional x-yplane, which represents the actual geographic wiring diagram as thefoundation of the visualization. The node inertia and the information ofnode geographic location are, obtained by fusion calculation on thetwo-dimensional plane to realize the visualization of the system inertiaspatial distribution. The inertia of each synchronous unit in thisscenario is shown in Table 3, which is based on the rated capacity ofthe generator itself, with a reference frequency f_(N)=50 Hz. Theestimation result of each node inertia in the grid is shown in Table 4,and the visual display diagram is shown in FIG. 4.

TABLE 3 Parameters of the generator models Generator H x′_(d) x′_(q)x_(d) x_(q) T′_(d0) T′_(q0) x_(l) 1 4.6 0.060 0.080 0.200 0.19 7.00 0.700.030 2 3.5 0.697 1.700 2.950 2.82 6.56 1.50 0.350 3 3.7 0.531 0.8762.495 2.37 5.70 1.50 0.304 4 2.8 0.436 1.660 2.620 2.58 5.69 1.50 0.2955 2.6 1.320 1.660 6.700 6.20 5.40 0.44 0.540 6 2.5 0.500 0.814 2.5402.41 7.30 0.40 0.224 7 2.6 0.490 1.860 2.950 2.92 5.66 1.50 0.322 8 4.60.570 0.911 2.900 2.80 6.70 0.41 0.280 9 3.1 0.570 0.587 2.106 2.05 4.791.96 0.298 10 4.5 0.310 0.080 1.000 0.69 10.20 0 0.125

TABLE 4 Estimation results of node inertia in the grid Node Node numberinertia/s 1 4.18 2 3.84 3 3.39 4 3.29 5 3.35 6 3.34 7 3.39 8 3.41 9 3.9410 3.31 11 3.32 12 3.28 13 3.28 14 3.20 15 2.93 16 2.83 17 3.00 18 3.1319 2.65 20 2.70 21 2.70 22 2.59 23 2.62 24 2.80 25 3.82 26 3.19 27 3.0928 3.06 29 3.04 30 2.70 31 4.80 32 3.32 33 3.51 34 2.62 35 2.56 36 2.5637 2.50 38 4.53 39 2.92

As illustrated in FIG. 4, the inertia of the left region of the systemis greater than 4 s, the inertia of the right region is less than 3 s,and the inertia of the middle region is between 3 s-4 s. Thedistribution of the system inertia gradually decreases from the upperleft to the lower right, which is consistent with the distribution ofthe actual system inertia, verifying the effectiveness of method forestimating the spatial distribution of the inertia based onmulti-innovation identification.

As illustrated in FIG. 5, based on the visual display diagram of theinertia level of the regional power grid obtained by the conventionalmethods, the power system is divided into four regions, and theequivalent inertia of each region is solved as follows. The equivalentinertia of region 1 is 4.47 s. The equivalent inertia of region 2 is 3.5s. The equivalent inertia of region 3 is 3.6 s. The equivalent inertiaof region 4 is 2.63 s. Compared with FIGS. 4 and 5, it shows that takingthe region as a research object in the spatial distribution of theinertia merely reflects the average situation of equivalent inertia in acertain region, which is difficult to show the internal spatialdistribution in detail. However, in this disclosure, the method forestimating the spatial distribution of the inertia based onmulti-innovation identification estimates the inertia values of allnodes in the system to clarify the specific situation of the spatialdistribution of system inertia, providing more effective information forthe dispatch and operation personnel to dig out the regions with weakinertia.

In summary, in this disclosure, the method for estimating the spatialdistribution of the inertia based on multi-innovation identificationestimates the overall spatial distribution of the inertia of the powersystem by estimating the inertia value of each node in the system,thereby refining the analysis of the spatial distribution of the inertiato a node level. And during the process of estimating the spatialdistribution of the inertia of the system, this method is almost free ofthe impact by the disturbance intensity and disturbance location, and isof good adaptability to different disturbances and strong robustness.Meanwhile, the time that the method based on multi-innovationidentification takes to solve the model is milliseconds, which iscapable of meeting the requirement of the real-time analysis for actualoperation of the power system during the analytical process, andeffectively help the dispatching personnel accurately dig out theweak-inertia regions, facilitating improving the adaptability of theaccess of the regional renewable energy/DC and the operation safety andstability of the system.

Embodiment 2

As illustrated in FIGS. 1-5, a device for implementing the method ofEmbodiment 1 is provided herein. The device includes an acquiring unit,a first calculating unit, a second calculating unit, and an analyzingunit. The acquiring unit is configured to acquire a time-series data andpreprocess the time-series data, in which the time-series data includesa frequency of each node of a power system and a transmission power of atransmission line. The acquiring unit is also configured to decompose afrequency signal on a series of different spaces to identify a momentwhen a disturbance occurs, and fast trace an impact of the disturbance.

The first calculating unit is configured to analyze a relationship amonga node inertia, a node frequency and a node injection power, constructan output error moving average (OEMA) model to describe anactive-power-frequency dynamic process of each node after thedisturbance, and analyze a coupling relationship between an unknownparameter of the OEMA model and the node inertia.

The second calculating unit is configured to solve the OEMA model in thefirst calculating unit by adopting a multi-innovation identification toidentify the unknown parameter of the OEMA model.

The analyzing unit is configured to calculate an inertia of each nodeusing the unknown parameter of the OEMA model identified by the secondcalculating unit according to an analytical relationship between avector of the unknown parameter and the node inertia to estimate anoverall spatial distribution of the inertia of the power system.

In this embodiment, the first calculating unit is performed as follows.

In terms of a mechanism of a frequency response of the power system, acoupling relationship among the node inertia, the node frequency and thenode injection power is analyzed.

Laplace transform is used to convert a time domain relationship amongthe node inertia, the node frequency, and the node injection power intoa continuous frequency domain expression, in which the continuousfrequency domain expression is a continuous transfer function includingan information of the node inertia.

The OEMA model is constructed to describe the active power-frequencydynamic change process of each node. Solving of the node inertia isconverted into identification of the unknown parameter of the OEMAmodel. And the coupling relationship between the unknown parameter ofthe OEMA model and the node inertia is analyzed. An expression of theOEMA model is shown as follows:

$\begin{matrix}{{y(t)} = {{\frac{B(z)}{A(z)}{u(t)}} + {{D(z)}{v(t)}}}} \\{{{x(t)} = {\frac{B(z)}{A(z)}{u(t)}}};}\end{matrix}$

where A(z)=1+a₁z⁻¹+ . . . +a_(n) _(a) z^(−n) ^(a) ; B(z)=b₁z⁻¹+ . . .+b_(n) _(b) z^(−n) ^(b) ; D(z)=1+d₁z⁻¹+ . . . +d_(n) _(d) z^(−n) ^(d) ;y(t) is an output sequence of the OEMA model to be identified; u(t) isan input sequence of the OEMA model to be identified; v(t) is azero-mean Gaussian random white noise sequence; x(t) is an internalvariable, which represents an actual undetectable variable of the powersystem; z⁻¹ is a unit backshift operator, expressed as z⁻¹y(t)=y(t−1);and a_(n) _(a) , b_(n) _(b) , d_(n) _(d) and represent parameters to beidentified of the OEMA model.

The estimation of the spatial distribution of the power system inertiabased on multi-innovation identification is performed in accordance withthe method in Embodiment 1, which will not be repeated herein.

Compared with the conventional devices for estimating the power systeminertia, the device of this disclosure refines the estimation result ofthe spatial distribution of the inertia to the node level, which has ahigher degree of refinement. This disclosure is capable of effectivelyovercoming the interference of a large number of random factors such asdata measurement noise, allowing for a strong robustness of the model.The multi-innovation identification adopted in this disclosure iscapable of effectively solving the model for estimating the spatialdistribution of the inertia. In this disclosure, the model is ofdesirable adaptability to disturbance faults. The model is capable ofaccurately estimating the spatial distribution of the inertia of thepower system under disturbances varying in location and intensity).

The above embodiments are only used to illustrate the objectives,technical solutions, and beneficial effects of this disclosure, and arenot intended to limit the present disclosure. It should be understoodthat any modifications, replacements and improvements made by thoseskilled in the art without departing from the spirit of this applicationshould fall within the scope of this application defined by the appendedclaims.

What is claimed is:
 1. A method for estimating an inertia distributionof a power system, the method comprising: (S1) acquiring time-seriesdata which comprises a frequency at each of a plurality of nodes of thepower system, and a transmission power of a transmission line, andpreprocessing the time-series data; decomposing the frequency in aseries of different frequency domains to identify a moment when adisturbance occurs and trace an impact of the disturbance fast; (S2)constructing an output error moving average (OEMA) model to describe anactive power-frequency dynamic process for each of the nodes of thepower system after the disturbance, and analyzing a couplingrelationship between an unknown parameter in the OEMA model and aninertia at each of the nodes of the power system; (S3) solving the OEMAmodel by adopting a multi-innovation identification, and identifying theunknown parameter of the OEMA model; and (S4) calculating the inertia ateach of the nodes of the power system by using the unknown parameter ofthe OEMA model identified in step (S3) according to an analyticalrelationship between a vector of the unknown parameter and the inertiaat each of the nodes of the power system, to estimate the inertiadistribution of the power system; wherein the step (S1) comprises: (S11)receiving, from a plurality of phasor measurement units (PMUs), thetime-series data; (S12) preprocessing the time-series data obtained instep (S11) through interpolation of missing data and elimination of adata measurement noise via-filtering; wherein the step (S2) comprises:(S21) in terms of a mechanism of a frequency response of the powersystem, analyzing a coupling relationship among the inertia at each ofthe nodes of the power system, the frequency at each of the nodes of thepower system and an injection power at each of the nodes of the powersystem; (S22) using Laplace transform to convert a time domainrelationship among the inertia, the frequency and the injection power ateach of the nodes of the power system into a continuous frequency domainexpression that is a continuous transfer function; and (S23)constructing the OEMA model to describe the active power-frequencydynamic change process for an identification of the unknown parameter ofthe OEMA model; and analyzing the coupling relationship between theunknown parameter of the OEMA model and the inertia at each of the nodesof the power system; wherein an expression of the OEMA model is shown asfollows: $\begin{matrix}{{y(t)} = {{\frac{B(z)}{A(z)}{u(t)}} + {{D(z)}{v(t)}}}} \\{{{x(t)} = {\frac{B(z)}{A(z)}{u(t)}}};}\end{matrix}$ wherein, A(z)=1+a₁z⁻¹+ . . . +a_(n) _(a) z^(−n) ^(a) ,B(z)=b₁z⁻¹+ . . . +b_(n) _(b) z^(−n) ^(b) , D(z)=1+d₁z⁻¹+ . . . +d_(n)_(d) z^(−n) ^(d) ; y(t) is an output sequence; u(t) is an inputsequence; v(t) is a zero-mean Gaussian random white noise sequence, x(t)is an internal variable, which represents an actual undetectablevariable of the power system; z⁻¹ is a unit backshift operator expressedas z⁻¹y(t)=y(t−1); and a_(n) _(a) , b_(n) _(b) , d_(n) _(d) representsparameters to be identified; wherein the step (S3) comprises: (S31)expanding a single innovation describing the active power-frequencydynamic change process into multiple innovations to establish amulti-innovation matrix; wherein the multi-innovation matrix isexpressed as follows: ${\begin{bmatrix}{X\left( {p,t} \right)} \\{Y\left( {p,t} \right)} \\{V\left( {p,t} \right)} \\{\psi\left( {p,t} \right)}\end{bmatrix} = \begin{bmatrix}{{x(t)},{x\left( {t - 1} \right)},\ldots,{x\left( {t - p + 1} \right)}} \\{{y(t)},{y\left( {t - 1} \right)},\ldots,{y\left( {t - p + 1} \right)}} \\{{v(t)},{v\left( {t - 1} \right)},\ldots,{v\left( {t - p + 1} \right)}} \\{{\varphi(t)},{\varphi\left( {t - 1} \right)},\ldots,{\varphi\left( {t - p + 1} \right)}}\end{bmatrix}};$ wherein φ(t) is an information vector comprising theinternal variable x(t), and φ(t)=[−x(t−1), . . . , −x(t−n_(a)), u(t−1),. . . , u(t−n_(b)), v(t−1), . . . , v(t−n_(d))]^(T); X(p,t) is anaccumulated internal variable; Y(p,t) is an accumulated output vector;V(p,t) is an accumulated noise variable; ψ(p,t) is an accumulatedinformation matrix; and p is the number of an innovation vector; (S32)replacing the internal undetectable variable of the OEMA model with anoutput of an auxiliary model; (S33) based on step (S32), setting acriterion function, obtaining an optimal identification result of aparameter of the OEMA model on a basis of a least squares estimationprinciple on the premise of allowing the criterion function to beminimum to obtain values of all unknown parameters of the OEMA model;wherein the step (S4) comprises: (S41) establishing a mappingrelationship between a z-plane and an s-plane based on bilineartransformation; and according to the vector of the parameter of the OEMAmodel identified in step (S3), obtaining an undetermined parameter of acontinuous transfer function corresponding to the OEMA model; (S42)performing an inverse Laplace transform on the transfer functionobtained in step (S41) to obtain an impulse response of a correspondingcontrol system; taking an initial response value to obtain a reciprocalof twice the inertia at each of the nodes of the power system, therebyestimating a value of the inertia at each of the nodes of the powersystem; (S43) visually displaying the inertia distribution of the powersystem by combining geographic information of each of the nodes, basedon the value of the inertia at each of the nodes of the power system,and locating a node that has lower inertia than other nodes of the powersystem.
 2. The method of claim 1, wherein the step (S1) comprises: (S13)normalizing the time-series data preprocessed in step (S12) to form adataset in a standard form; (S14) decomposing the frequency at each ofthe nodes of the power system in the series of different frequencydomains based on wavelet multi-resolution analysis using a DaubechiesWavelet; and extracting a wavelet detail coefficient at each level todetermine a modulus maxima spot, thereby detecting a singular point ofthe frequency at each of the nodes of the power system that is aninitial time when the disturbance occurs.
 3. The method of claim 1,wherein in step (S21), the coupling relationship among the node inertia,the node frequency and the node injection power is expressed as follows:${{2H_{i}\frac{d\Delta{f_{i}(t)}}{dt}} = {\Delta{P_{i}(t)}}};$ which isfurther converted into:${H_{i} = \frac{\Delta{P_{i}(t)}}{2\frac{d\Delta{f_{i}(t)}}{dt}}};$wherein H_(i) is the value of the inertia at a node i, and represents amagnitude of a hindering effect of the inertia on the frequency at thenode i during a ΔP_(i)−dΔfi/dt process; Δf_(i) is a frequency deviation;ΔP_(i) is an increment of the injection power; in step (S22), thecontinuous transfer function corresponding to the relationship among theinertia at each of the nodes of the power system, the frequency at eachof the nodes of the power system and the injection power at each of thenodes of the power system is expressed as follows:${{G(s)} = {\frac{\Delta{f_{i}(s)}}{\Delta{P_{i}(s)}} = \frac{1}{2H_{i}s}}};$wherein Δf_(i)(s) is a frequency domain expression of the frequencydeviation; ΔP_(i)(s) is a frequency domain expression of the incrementof the injection power at each of the nodes of the power system; and sis a Laplace operator.
 4. The method of claim 1, wherein in step (S41),the bilinear transform is configured to realize a mapping from thez-plane to the s-plane, and is expressed as follows:${z = {\frac{e^{s{T/2}}}{e^{{- s}{T/2}}} \approx \frac{1 + {s{T/2}}}{1 - {s{T/2}}}}};$wherein z is a sampling Laplace operator; s is the Laplace operator; andTis a sampling period.
 5. The method of claim 1, wherein the frequencyat each of the nodes of the power system is decomposed and reconstructedunder the Daubechies Wavelet basis to detect the singular point of thefrequency to identify an initial moment upon a response of the inertiaat each of the nodes of the power system, wherein an expression of thewavelet multi-resolution analysis is shown as follows:${{W_{f}\left( {a,b} \right)} = {< f}},{{\psi_{a,b}>={{❘a❘}^{{- 1}/2}{\int_{- \infty}^{+ \infty}{{f(t)}{\psi\left( \frac{t - b}{a} \right)}{dt}}}}};}$wherein ψa, b(t) is a mother wavelet function, which is the DaubechiesWavelet basis; a is a scale coefficient; b is a translating coefficient,and f(t) is an acquired frequency.